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Abstract 

We present theoretical studies of the influence of spin orbit coupling on the spin wave excitations 
of the Fe monolayer and bilayer on the W(110) surface. The Dzyaloshinskii-Moriya interaction is 
active in such films, by virtue of the absence of reflection symmetry in the plane of the film. When 
the magnetization is in plane, this leads to a linear term in the spin wave dispersion relation for 
propagation across the magnetization. The dispersion relation thus assumes a form similar to 
that of an energy band of an electron trapped on a semiconductor surfaces with Rashba coupling 
active. We also show SPEELS response functions that illustrate the role of spin orbit coupling in 
such measurements. In addition to the modifications of the dispersion relations for spin waves, the 
presence of spin orbit coupling in the W substrate leads to a substantial increase in the linewidth 
of the spin wave modes. The formalism we have developed applies to a wide range of systems, 
and the particular system explored in the numerical calculations provides us with an illustration 
of phenomena which will be present in other ultrathin ferromagnet/substrate combinations. 
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I. INTRODUCTION 



The study of spin dynamics in ultrathin ferromagnets is of fundamental interest, since new 
physics arises in these materials that has no counterpart in bulk magnetism. Examples are 
provided by relaxation mechanisms evident in ferromagnetic resonance and Brillouin light 
scattering studies,- - - and also for the large wave vectors probed by spin polarized electron 
loss spectroscopy (SPEELS).- Of course, by now the remarkable impact of ultrathin film 
structures on magnetic data storage is very well known, and other applications that exploit 
spin dynamics in such materials are envisioned. Thus these issues are important from a 
practical point of view as well as from that of fundamental physics. 

Theoretical studies of the nature of spin waves in ultrathin films adsorbed on metal sub- 
strates have been carried out for some years now, along with comparison with descriptions 
provided with the Heisenberg model.- In this paper, we extend the earlier theoretical treat- 
ments to include the influence of spin orbit coupling on the spin wave spectrum of ultrathin 
films. This extension is motivated by a most interesting discussion of the ground state of the 
Mn monolayer on the W(110) surface. A nonrelativistic theoretical study of this system pre- 
dicted that the ground state would be antiferromagnetic in character.- This prediction was 
confirmed by spin polarized scanning tunneling microscope studies of the system.- However, 
recent experimental STM data with a more sensitive instrument showed a more complex 
ground state, wherein the ground state is in fact a spin density wave.- One can construct 
the new state by beginning with the antiferromagnet, and then superimposing on this a 
long wavelength modulation on the direction of the moments on the lattice. The authors 
of ref. |9| argued that the lack of reflection symmetry of the system in the plane of the film 
activates the Dzyaloshinskii Moriya (DM) interaction, and the new state has its origin in this 
interaction. They also presented relativistic and ab initio calculations that gave an excellent 
account of the new data. The reflection symmetry is broken simply by the presence of the 
substrate upon which the film is grown. This argument to us is most intriguing, since one 
can then conclude that the DM interaction must be active in any ultrathin ferromagnet; 
the substrate is surely always present. The DM interaction has its origin in the spin orbit 
interaction, which of course is generally very weak in magnets that incorporate the 3d tran- 
sition elements as the moment bearing entities. However, in the case of the Mn monolayer 
on W(110) hybridization between the Mn 3d and the W 5d orbitals activates the very large 
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W spin orbit coupling, with the consequence that the strength of the DM interaction can 
be substantial, as illustrated by the calculations presented in ref. Ill One may expect to see 
substantial impact of the DM interaction in other ultrathin magnets grown on 5d substrates, 
and possibly 4d substrates as well. 

We have here another example of new physics present in ultrathin magnets that is not 
encountered in the bulk form of the material from which the ultrathin structure is fabricated. 
The purpose of this paper is to present our theoretical studies of spin orbit effects on spin 
waves and also on the dynamic susceptibility of a much studied ultrathin film/substrate 
combination, the Fe monolayer and bilayer on W(110). We find striking effects. For instance, 
when the magnetization is in plane, as we shall see the DM interaction introduces a term 
linear in wave vector in the dispersion relation of spin waves. Thus the uniform spin wave 
mode at zero wave vector acquires a finite group velocity. We find this to be in the range 
of 2 x 10 5 cm/sec for the Fe monolayer on W(110). Furthermore, left/right asymmetries 
appear in the SPEELS response functions. Thus, we shall see that spin orbit coupling has 
clear effects on the spin excitations of transition metal ultrathin ferromagnets grown on 5d 
substrates. 

We comment briefly on the philosophy of the approach used here, and in various earlier 
publications.- Numerous authors proceed as follows. One may generate a description of the 
magnetic ground state of the adsorbed films by means of an electronic structure calculation 
based on density functional theory. It is then possible to calculate, within the framework of 
an adiabatic approximation, effective Heisenberg exchange integrals Jij between the magnetic 
moments in unit cell i and unit cell j. These may be entered into a Heisenberg Hamilto- 
nian, and then spin wave dispersion relations may be calculated through use of spin wave 
theory. It has been known for decades^ that in the itinerant 3d magnets, effective exchange 
interactions calculated in such a manner have very long range in real space. Thus, one must 
include a very large number of distant neighbors in order to obtain converged results. This 
is very demanding to do with high accuracy for the very numerous distant neighbors, since 
the exchange interactions become very small as one moves out into distant neighbor shells. 

At a more fundamental level, as noted briefly above, discussions in earlier publications 
show that in systems such as we study here, the adiabatic approximation breaks down badly, 
with qualitative consequences.- First, spin wave modes of finite wave vector have very short 
lifetimes, by virtue of decay into the continuum of particle hole pairs (Stoner excitations) 
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even at the absolute zero of temperature*^ whereas in Heisenberg model descriptions their 
lifetime is infinite. In multi layer films, the earlier calculations show that as a consequence of 
the short lifetime, the spectrum of spin fluctuations at large wave vectors contains a single 
broad feature which disperses with wave vector in a manner similar to that of a spin wave; 
this is consistent with SPEELS data on an eight layer film of Co on Cu(100). 4 This picture 
stands in contrast to that offered by the Heisenberg model, in which a film of N layers has 
N spin wave modes for each wave vector, and each mode has infinite lifetime. 

The method developed earlier, and extended here to incorporate spin orbit coupling, 
takes due account of the breakdown of the adiabatic approximation and also circumvents 
the need to calculate effective exchange interactions in real space out to distant neighbor 
shells. We work directly in wave vector space through study of the wave vector and frequency 
dependent susceptibility discussed below, denoted as x+,-(Q\\> l±, l'±)- The imaginary part 
of this object, evaluated for l± = l' ± and considered as a function of frequency Q for fixed 
wave vector Q\\ provides us with the frequency spectrum of spin fluctuations on layer l± for 
the wave vector chosen. Spin waves appear as peaks in this function, very much as they do 
in SPEELS data, and in a manner very similar to that used by experimentalist we extract 
a dispersion relation for spin waves by following the wave vector dependence of the peak 
frequency. We never need to resort to a real space summation procedure over large number 
of neighbors, coupled by very tiny exchange couplings. The spin wave exchange stiffness 
can be extracted either by fitting the small wave vector limit of the dispersion relation so 
determined, or alternatively by utilizing an expression derived earlier- which once again does 
not require a summation in real space. 

We comment on another feature of the present study. In earlier calculations,- 1 ^^ as in the 
present paper, an empirical tight binding description forms the basis for our description of 
the electronic structure. Within this approach, referred to as a multi band Hubbard model, 
we can generate the wave vector and frequency dependent susceptibility for large systems. 
In the earlier papers, effective tight binding parameters were extracted from bulk electronic 
structure calculations. The present studies are based on tight binding parameters obtained 
directly from a RS-LMTO-ASA calculation for the Fe/W(110) system. We also obtain tight 
binding parameters by fitting KKR based electronic structure calculations for the ultrathin 
film/substrate combinations of interest. We find that spin waves in the Fe/W(110) system 
are quite sensitive to the empirical tight binding parameters which are employed, though as 



4 



we shall see the various descriptions provide very similar pictures of the one electron local 
density of states. 

We note that Udvardi and Szunyogh 12 have also discussed the influence of spin orbit 
coupling on the dispersion relation of spin waves in the Fe monolayer on W(110) within 
the framework of the adiabatic approach discussed above, where exchange interactions and 
other magnetic parameters are calculated in real space. We shall discuss a comparison with 
our results and theirs below. There are differences. Most particularly, we note that in 



Fig. 3, the authors of ref. [12| provide two dispersion curves for propagation perpendicular 
to the magnetization, whereas in a film such as this with one spin per unit cell there can 
be only one magnon branch. Additionally and very recently, Bergmann and coworkers^ 
investigated within an adiabatic approach finite temperature effects on the magnon spectrum 
of Fe/W(110). 

In section II, we comment on our means of introducing spin orbit coupling into the theory. 
The results of our calculations are summarized in Section III and concluding remarks are 
found in section IV. 



II. CALCULATION OF THE DYNAMIC SUSCEPTIBILITY IN THE PRESENCE 
OF SPIN ORBIT COUPLING 

The formalism for including spin orbit coupling effects in our description of spin dynamics 
is quite involved, so in this section we confine our attention to an outline of the key steps, 
and an exposition of the overall structure of the theory. Our starting point is the multi band 
Hubbard model of the system that was employed in our earlier study of spin dynamics in 
ultrathin ferromagnets. The starting Hamiltonian is written as^ 

ij fil/a /ivfi'v' icrcr' 

where i and j are site indices, a, a' refer to spin, and fi, v to the tight binding orbitals, nine in 
number for each site, which are included in our treatment. The Coulomb interactions operate 
only within the 3d orbitals on a given lattice site. The film, within which ferromagnetism 
is driven by the Coulomb interactions, sits on a semi-infinite substrate within which the 
Coulomb interaction is ignored. 
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In our empirical tight binding picture, the spin orbit interaction adds a term we write as 
H so = J2 J2 7T [ L ^( c lt c ^t - 4^ iH ) + L^cj^Ci^ + L^c\^c iH ] (2) 

where L is the angular momentum operator, Aj is the local spin-orbit coupling constant, 
L ± = L x ±iL y and L° u = (n\L a \i>). We assume that the spin orbit interaction, present both 
within the ferromagnetic film and the substrate, operates only within the 3d atomic orbitals. 
A convenient tabulation of matrix elements of the orbital angular momentum operators is 
found in ref. 13. 



Information on the spin waves follows from the study of the spectral density of the trans- 
verse dynamic susceptibility x+-(Q\\M'J±J'±) as discussed above. From the text around 



Eq. (1) of ref. [14], we see that this function describes the amplitude of the transverse spin 
motion (the expectation value of the spin operator S + in the layer labeled l±) to a fictitious 
transverse magnetic field of frequency Q and wave vector Q\\ parallel to the film surface that 
is applied to layer l' ± of the sample. The spectral density, given by lm{x+-(Q\\ , Q', l±, 
when multiplied by the Bose Einstein function n(Q) = [exp(/3f2) — l]^ 1 is also the amplitude 
of thermal spin fluctuations of wave vector Q\\ and frequency f2 in layer l±. We obtain in- 
formation regarding the character (frequency, linewidth, and amplitude in layer l±) of spin 
waves from the study of this function, as discussed earlier.- 

Our previous analyses are based on the study of the dynamic susceptibility just described 
through use of the random phase approximation (RPA) of many body theory. The Feynman 
diagrams included in this method are the same as those incorporated into time dependent 
density functional theory, though use of our Hubbard model allows us to solve the result- 
ing equation easily once the very large array of irreducible particle hole propagators are 
generated numerically. 

Our task in the present paper is to extend the RPA treatment to incorporate spin orbit 
coupling. The extension is non trivial. The quantity of interest, referred to in abbreviated 
notation as x+-i ma y be expressed as a commutator of the spin operators S + and S~ whose 
precise definition is given earlier.-*^ With spin orbit coupling ignored, the RPA decoupling 
procedure leads to a closed equation for x+,-- When the RPA decoupling is carried out in 
its presence, we are led to a sequence of four coupled equations which include new objects 
we may refer to as X-,-i Xt- an d Xi-- The number of irreducible particle hole propagators 
that must be computed likewise is increased by a factor of four. For a very simple version 



of a one band Hubbard model, and for a very different purpose, Fulde and Luther carried 
out an equivalent procedure many years ago^. In what follows, we provide a summary of 
key steps along with expressions for the final set of equations. 

To generate the equation of motion, we need the commutator of the operator S^ v (l, I') = 
c Lt c i'j4 with the Hamiltonian. One finds 

[S+M, l'),H S o] = \ XX A ' ,L i c lt c ^t - X i L w c li c i'H + \>L z uv c U c i'vl ~ X i LZ r,A c i'^}- ( 3 ) 
z v 

The last two terms on the right hand side of Eq. [3] lead to terms in the equation of motion 
which involve x+- whereas the first two terms couple us to the entities Xf,- an d Xl,-- 
When we write down the commutator of these new correlation functions with the spin orbit 
Hamiltonian, we are led to terms which couple into the function X-,- which is formed from 
the commutator of two S~~ operators. In the absence of spin orbit coupling, a consequence of 
spin rotation invariance of the Hamiltonian is that the three new functions just encountered 
vanish. But they do not in its presence, and they must be incorporated into the analysis. 

One then introduces the influence of the Coulomb interaction into the equation of motion, 
and carries out an RPA decoupling of the resulting terms. The analysis is very lengthy, so 
here we summarize only the structure that results from this procedure. Definitions of the 
various quantities that enter are given in the Appendix. We express the equations of motion 
in terms of a 4 x 4 matrix structure, where in schematic notation we let x = X+,—> 
X — Xt,-> X^ — Xi- an d X^ — X--- The four coupled equations then have the form 

q x (s) = A (s) + J2(B SS ' + B ss ') X {s,) (4) 

s' 

Each quantity in Eq. H] has attached to it four orbital indices, and four site indices. To be 
explicit, x^ = Xt>- which enters Eq. (4) is formed from the commutator of the operator 
c Lt c ' /l/ t with c^/^CmV't an d in full we denote this quantity as X^um'v'i^''' mm ')- The site 
indices label the planes in the film, and we suppress reference to Q and Q\\. The products 
on the right hand side of Eq. H] are matrix multiplications that involve these various indices. 
For instance, the object J2s> B ss x^ is labeled by four orbital and four site indices so 

[B ss 'x {s,) UMtt'i mm') B^ s (ll>; nn') X %Ann'; mm'). (5) 

7<5 nn' 

One proceeds by writing Eq. H] in terms of the dynamic susceptibilities that characterize 
the non-interacting system. These, referred to also as the irreducible particle hole propaga- 
tors, are generated by evaluating the commutators which enter into the definition of x^ in 



the non interacting ground state. These objects, denoted by obey a structure similar 
to Eq. HI 

Qx (0s) = A ^ + J2B ss 'x {s '\ (6) 

s> 

It is then possible to relate x to x^ 0s ' 1 through the relation, using four vector notation, 

m = x (0) m + (n-B)- i Bm- (7) 

The matrix structure T = (Q — B)^ 1 may be generated from the definition of B, which may 
be obtained from the equation of motion of the non-interacting susceptibility, Eq. [6j Then 
B follows from the equation of motion of the full susceptibility, as generated in the RPA. 
One may solve Eq. [7] 

x(n) = [i-(n-B)- 1 B]- 1 x {0 \n), (8) 

so our basic task is to compute the non interacting susceptibility matrix x^ and then carry 
out the matrix inversion operation displayed in Eq. [HJ For this we require the single particle 
Greens functions (SPGFs) associated with our approach. 

To generate the SPGFs, we set up an effective single particle Hamiltonian H sp by intro- 
ducing a mean field approximation for the Coulomb interaction. The general structure of 
the single particle Hamiltonian is 

Hsp — ^ ^ T^j Cifj^Cjva + 5Z{ a i;/ji/ C i/4 C ^t + a i;^ C lfi^ C ivl} (9) 
ij fJ-fcr i flu 

where the effective hopping integral T^ ucr contains the spin diagonal portion of the spin orbit 
interaction, along with the mean field contributions from the Coulomb interaction. The form 
we use for the latter is stated below. The coefficients in the spin flip terms are given by 

OL%iw = — ^2 Ui- V n tlll (c iv ^Ciyf) . (10) 

We then have the eigenvalue equation that generates the single particle eigenvalues and 
eigenfunctions in the form H sp \<p s ) = E s \(j) s )\ we can write this in the explicit form 

J2 [Saa'Tu 1 "' + duiS^^a*.^ + S^S^auw)] (VVs) = E s (i/ia\<j) s }. (11) 

l r\a' 

The single particle Greens function may be expressed in terms of the quantities that enter 
Eq. [HJ We have for this object the definition 

GW;,w(£) = -^(t)({c^(t), 4^,(0)}) (12) 
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and one has the representation 

~ \ n-E s + iV ■ (13) 

These functions may be constructed directly from their equations of motion, which read, 
after Fourier transforming with respect to time, 



l ncr" 

(14) 

For the case where the substrate is semi infinite, our means of generating a numerical solution 
to the hierarchy of equations stated in Eq. [TH has been discussed earlier. What remains is 
to describe how the Coulomb interaction enters the effective hopping integrals T^j ucr that 
appear in Eq. [9j Eq. [11] and Eq. [TH 

There are, of course, a large number of Coulomb matrix elements in the original Hamil- 
tonian, even if the Coulomb interactions are confined to within the 3d shell. Through the 
use of group theory— the complete set of Coulomb matrix elements may be expressed in 
terms of three parameters. These are given in Table I of the first cited paper in ref. 0. 
In subsequent work, we have found that a much simpler structure^ nicely reproduces re- 
sults obtained with the full three parameter form. We use the simpler one parameter form 
here, for which E/j iAM/iM v = Ui5^ u >5^ u . Then in the mean field approximation, the Coulomb 
contribution to the single particle Hamiltonian assumes the form 

h£> = -Y^r Bit^t - (is) 

Here is the magnitude of the moment on site i. The Coulomb interactions U,i are non zero 
only within the ultrathin ferromagnet, and the moments rrii, determined self consistently, 
vary from layer to layer when we consider multi layer ferromagnetic films. 

It should be noted that when the Ansatz just described is employed in Eq. [TO], the term 
from the Coulomb interaction on the right hand side becomes proportional to the transverse 
component of the moment located on site i and this vanishes identically. Thus, despite the 
complexity introduced by the spin orbit coupling, when the simple one parameter Ansatz for 
the Coulomb matrix elements is employed, one needs no parameters beyond the moment on 
each layer in the self consistent loops that describe the ground state. In the present context, 
this is an extraordinarily large savings in computational labor, and this will allow us to 
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address very large systems in the future. It is the case that certain off diagonal elements 
such as {(? m u'\. c i'ir\) appear in the quantities defined in the Appendix. Notice, for example, 
the expressions in Eqs. IA.1I Once the ground state single particle Greens functions are 
determined, such expectation values are readily computed. 

III. RESULTS AND DISCUSSION 

In earlier studies of Fe layers on W(110),— ^ as noted above, the electronic structure 
was generated through use of tight binding parameters obtained from bulk electronic struc- 
ture calculations. These calculations generate effective exchange interactions comparable 
in magnitude to those found in the bulk transition metals,— with the consequence that for 
both monolayer Fe and bilayer Fe on W(110) the large wave vector spin waves generated 
by theory are very much stiffer than found experimentally 2 ^* 2 ^ though it should be noted 
that for the bilayer, the calculated value of the spin wave exchange stiffness is in excellent 
accord with the data. 23 Subsequent calculations which construct the spin wave dispersion 
relation from adiabatic theory based on calculations of effective exchange integrals also gen- 
erate spin waves for the monolayer substantially stiffer than found experimentally,— though 
they are softer than in our earlier work by a factor of two or so. We remark that it has 
been suggested that the remarkably soft spin waves found experimentally may have origin 
in carbon contamination of the monolayer and bilayer.— We remark here that this can be 
introduced during the SPEELS measurement. We note that the magnetic properties of Fe 
monolayers grown on carbon free W(110) 24 differ dramatically from those grown on surfaces 
now known to be contaminated by carbon.— In the former case, the domain walls have a 
thickness of 2.15 nm,— whereas in the latter circumstance very narrow walls with thickness 
bounded from above by 0.6 nm are found.™ This suggests that the strength of the effective 
exchange is very different in the two cases, with stiffer exchange in the carbon free samples. 

The considerations of the previous paragraph have motivated us to carry out a series of 
studies of the effective exchange in the Fe monolayer on W(110) within the framework of 
three different electronic structure calculations. We find that although all three give local 
density of states that are very similar, along with very similar energy bands when these are 
examined, the intersite exchange interactions vary substantially. First, we have employed the 
parameter set used earlier that is based on bulk electronic structures' 2 ^ in new calculations 
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we call case A. In case B, we have employed an approach very similar to that used in ref. 
though in what follows our calculation of effective exchange integrals is non relativistic. 
This is the Korringa Kohn Rostoker Greens Function (KKR-GF) method,— which employs 
the atomic sphere approximation and makes use of the Dyson equation G = g + gVG as 
given in matrix notation. This allows us to calculate the Greens function G of an arbitrary 
complex system given the perturbing potential V and the Greens function g of a reference 
unperturbed system. Within the Local Spin Density Approximation (LSDA), 27 We consider 
a slab of five monolayers of W with the experimental lattice constant on top of which an Fe 
monolayer is deposited and relaxed by -12. 9% 12 with respect to the W interlayer distance. 
Angular momenta up to / max = 3 were included in the Greens functions with a k mesh of 
6400 points in the full two dimensional Brillouin zone. The effective exchange interactions 
were calculated within the approximation of infinitesmal rotations 2 ^ that allows one to use 
the magnetic force theorem. This states that the energy change due to infinitesmal rotations 
in the moment directions can be calculated through the Kohn Sham eigenvalues. 

Method C is the Real Space Linear-Muffin-Tin-Orbital approach as implemented, also, 
in the atomic sphere approximation (RS-LMTO-ASA).— ~— Due to its linear scaling, this 
method allows one to address the electronic structure of systems with a large number of 
atoms for which the basic eigenvalue problem is solved in real space using the Haydock 
recursion method. The Fe overlayer on the W(110) substrate was simulated by a large 
bcc slab which contained ~6800 atoms, arranged in 12 atomic planes parallel to the (110) 
surface, with the experimental lattice parameter of bulk W. One empty sphere overlayer 
is included, and self consistent potential parameters were obtained for the empty sphere 
overlayer, the Fe monolayer, and the three W layers underneath using LSDA.— For deeper 
W layers we use bulk potential. Nine orbitals per site (the five 3d and 4 sp complex) were 
used to describe the Fe valence band and the empty sphere overlayer, and for W the fully 
occupied 4f orbitals were also included in the core. To evaluate the orbital moments we use 
a scalar relativistic (SR) approach and include a spin orbit coupling term XL ■ S at each 
variational step.— In the recursion method the continued fraction has been terminated after 
30 recursion levels with the Beer Pettifor terminator.— The TB parameters so obtained are 
inserted into our semi-empirical scheme and this allows us to generate the non interacting 
susceptibilities which enter our full RPA description of the response of the structure. 

In order to compare the electronic structures generated by the approaches just described, 
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we turn our attention to the local density of states for the majority and minority spins in 
the adsorbed Fe monolayer. These are summarized in Fig. [TJ 

The local densities of states (LDOS) generated by the three sets of TB parameters have 
approximately the same overall features, as we see from Fig. [TJ The main differences appear 
in the majority spin band, which overlaps the 5d states in the W substrate over a larger 
energy range than the minority band. This is also true if we compare the LDOS generated 
by the tight binding parameters extracted from the KKR electronic structure to the LDOS 
obtained directly from the KKR calculations (red dashed line in Fig.[Tb). The Fe-W hopping 
parameters are indeed the least accurate portion of our parametrization scheme. In case A 
we just used the Fe-Fe bulk parameters to describe the Fe-W hopping. In case B we extracted 
TB parameters for Fe by fitting a KKR calculation of an unsupported Fe monolayer with a 
lattice parameter matching that of the W substrate. For the Fe-W hopping we used the Fe 
parameters obtained from the fitting, scaled to mimic the Fe-W distance relaxation. The 
relaxation parameter was chosen to give the correct spin magnetic moment for the adsorbed 
Fe monolayer. In case C all parameters were directly provided by the RS-LMTO-ASA code, 
but in the DFT calculations the Fe-W distance was assumed to be equal to the distance 
between W layers. Thus, the main difference between cases B and C is the treatment of the 
mixing between Fe and W states and this is expected to affect more strongly the states that 
occupy the same energy range. 

As noted above, while the local density of states provided by the three approaches to 
the electronic structure are quite similar as we see from Fig. [TJ (and the same is true of the 
electronic energy bands themselves if these are examined), the exchange interactions differ 
substantially for the three cases. For the first, second and third neighbors we have (in meV) 
42.5, 3.72 and 0.46 for model A, 28.7, -7.87 and 0.31, for model B and 11.23, -7.31 and 0.22 
for model C. The authors of ref. [l2| find 10.84, -3.34 and 3.64 for these exchange integrals. 

We now turn to our studies of spin excitations in the Fe monolayer and Fe bilayer on 
W(110) within the framework of the electronic structure generated through use of the ap- 
proach in case C. We will discuss the influence of spin orbit coupling on both the trans- 
verse wave vector dependent susceptibility though study of the spectral density function 
A(Q\\,Q,;l±) = — lm.{x+-(Q\\, ^; l±)} discussed in section I. This function, for fixed wave 
vector Qn, when considered as a function of frequency Q, describes the frequency spectrum 
of the fluctuations of wave vector Qn of the transverse magnetic moment in layer l± as noted 
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FIG. 1: (color online) For the Fe monolayer adsorbed on W(110), we show the local density of 
states in the Fe monolayer. The majority spin density of states is shown positive and the minority 
spin density of states is negative. The zero of energy is at the Fermi energy. In (a), bulk electronic 
structure parameters are used as in the second of the two papers cited in |5| (caseA) . In (b) , we 
have the density of states generated by method B. The black curve is found by fitting the KKR 
electronic energy band structure to tight binding parameters as described in the text, and the red 
curve is calculated directly from the KKR calculation. In (c) we have the local density of states 
generated by method C. 
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above. In the frequency regime where spin waves are encountered, this function is closely 
related to (but not identical to) the response function probed in a SPEELS measurement. 

In Fig. [2J for the Fe monolayer on W(110), we show the spectral density function cal- 
culated for three values of \Q\\\, for propagation across the magnetization. Thus, the wave 
vector is directed along the short axis in the surface. This is the direction probed in SPEELS 
studies of the Fe monolayer on this surface.— In each figure, we show three curves. The green 
dashed curve is calculated with spin orbit coupling set to zero. We show only a single curve 
for this case, because the spectral density is identical for the two directions of propagation 
across the magnetization, +Q\\ and —Q\\- When spin orbit coupling is switched on, for the 
two directions just mentioned the response function is very different, as we see from the red 
and black curve in the various panels. These spin wave frequencies, deduced from the peak in 
the response functions as discussed in section I, differ for the two directions of propagation, 
and also note that the peak intensities and linewidths differ as well. It is the absence of both 
time reversal symmetry and reflection symmetry which renders +Q\\ and —Q\\ inequivalent 
for this direction of propagation. The system senses this breakdown of symmetry through 
the spin orbit interaction. If one considers propagation parallel to the magnetization, the 
asymmetries displayed in Fig. [2] are absent. The reason is that for this direction of propaga- 
tion, reflection in the plane that is perpendicular to both the magnetization and the surface 
is a good symmetry operation of the system, but takes +Q\\ into —Q\\ thus rendering the two 
directions equivalent. Recall, of course, that the magnetization is a pseudo vector in regard 
to reflections. Notice how very broad the curves are for large wave vectors; the lifetime of 
the spin waves is very short indeed. 

As discussed in section I, we may construct a spin wave dispersion curve by plotting the 
maxima in spectral density plots such as those illustrated in Fig. [2] as a function of wave 
vector. We show dispersion relations constructed in this manner in Fig. [3j with spin orbit 
coupling both present and absent. In Fig. [3^, and for propagation perpendicular to the 
magnetization we show the dispersion curve so obtained for wave vectors throughout the 
surface Brillouin zone, and in Fig. [3b we show its behavior for small wave vectors. 

Let is first consider Fig, 3(a). Here the dispersion curve extends throughout the two 
dimensional Brillouin zone. At the zone boundary, quite clearly the slope of the dispersion 
curve does not vanish. In this direction of propagation, the nature of the point at the zone 
boundary does not require the slope to vanish. What is most striking, clearly, is the anomaly 
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FIG. 2: (color online) The spectral density functions A{Qh , O; l±) evaluated in the Fe monolayer for 
three values of the wave vector in the direction perpendicular to the magnetization. We have (a) 
|Q||| = 0.4A , (b) |Q||| = l.OA and (c) |Q[|| = 1.4A" . The green curve (dashed) is calculated 
with spin orbit coupling set to zero; the spectral density here is independent of the sign of Qn. The 
red and black curves are calculated with spin orbit coupling turned on. Now we see asymmetries 
for propagation across the magnetization, with the red curve Qn directed from left to right and the 
black curve from right to left. 

in the vicinity of 1 A . This feature is evident in the calculation with spin orbit coupling 
absent, and for positive values of the wave vector the feature becomes much more dramatic 
when spin orbit coupling is switched on. Anomalies rather similar to those in the black 
curve in Fig. 3(a) appear in the green dispersion curve found in Fig. 3 of ref. 12, though 
these authors did not continue their calculation much beyond the 1 A 1 regime. Our s pin 



waves are very much softer than theirs in this spectral region, notice. In Fig. 3 of ref. 
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one finds two dispersion curves, one a mirror image of the second. Thus, these authors 
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FIG. 3: (color online) Spin wave dispersion relations constructed from peaks in the spectral den- 
sity, for the Fe monolayer on W(110). The wave vector is in the direction perpendicular to the 
magnetization. The red curve is constructed in the absence of spin orbit coupling, it is included in 
the black curve. 



display two spin wave frequencies for each wave vector. This surely is not correct. For a 
structure with one atom per unit cell, there is one and only one spin wave mode for each 
wave vector, though as discussed above for the structure explored here symmetry allow the 
left/right asymmetry in the dispersion curve illustrated in our Fig. [3j 

In Fig. |3b, again with spin orbit coupling switched on and off, we show an expanded 
view of the dispersion curve for small wave vectors. With spin orbit interaction switched off, 
at zero wave vector we see a zero frequency spin wave mode, as required by the Goldstone 
theorem when the underlying Hamiltonian is form invariant under spin rotation. The curve 
is also symmetrical, and is accurately fitted by the form Q(Q\\) = 149Q| (meV), with the 
wave vector in A \ whereas with spin orbit coupling turned on the dispersion relation is 
fitted by fi(Q\\) = 3.4 — 11.8Q|| + 143Qy (meV). Spin orbit coupling introduces an anisotropy 
gap at Q\\ = 0, and most striking is the term linear in wave vector. This has its origin in the 
Dzyaloshinskii Moriya interaction whose presence, as argued by the authors of ref. |9|, has its 
origin in the absence of both time reversal and inversion symmetry, for the adsorbed layer. 

At long wavelengths, one may describe spin waves by classical long wavelength phe- 
nomenology. The linear term in the dispersion curve has its origin in a term in the energy 
density of the spin system of the form 

dS„(x, z) 



V, 



DM 



dx 



(16) 
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Here S a (x, z) is a spin density, the xz plane is parallel to the surface, and the magnetization 
is parallel to the z direction. 

One interesting feature of the spin wave mode whose dispersion relation is illustrated in 
Fig. [3b is that at Q» = 0, the mode has a finite group velocity. The fit to the dispersion 
curve gives this group velocity to be —^q - ~ 2 x 10 5 cm/s , which is in the range of acoustic 
phonon group velocities. 

We turn now to our calculations of spin waves and the response functions for the Fe bilayer 
on W(110). Let us first note that experimentally the orientation of the magnetization in 
the bilayer appears to be dependent on the surface upon which the bilayer is grown. For 
instance, when the bilayer is on the stepped W(110) surface, it is magnetized perpendicular 
to the surface,— a result in agreement with ab initio calculations of the anisotropy realized in 
the epitaxial bilayer. 3 - However, in the SPEELS studies of spin excitations in the bilayer—^ 
the magnetization is in plane. In our calculations, we find for model B the magnetization is 
perpendicular to the surface, whereas in model C it lies in plane, along the long axis very 
much as in the SPEELS experiments. The anisotropy in the bilayer is not particularly large, 
on the order of 0.5 meV/Fe atom, and one sees from these results that it is a property quite 
sensitive to the details of the electronic structure. The fact that model B and model C give 
the two different stable orientation of the magnetization allows us to explore spin excitations 
for the two different orientations of the magnetization. 

We first turn our attention to the case where the magnetization lies in plane. The bilayer 
has two spin wave modes, an acoustic mode for which the magnetization in the two planes 
precesses in phase, and an optical mode for which they precess 180 degrees out of phase. 
In Fig. HI we show calculations of the dynamic susceptibility in the frequency range of the 
acoustic mode for two values of the wave vector, Q\\ = +0.5A 1 and Q\\ = — 0.5A 1. A 
spin orbit induced left right asymmetry is clearly evident both in the peak frequency and 
the height of the feature. Very recently, beautiful measurements of spin orbit asymmetries 
in the Fe bilayer have appeared,— and the results of our Fig. H] are to be compared with 
Fig. 3 of ref. |39|. Theory and experiment are very similar, both in regard to the intensity 
asymmetry and also the spin orbit induced frequency shift, though our calculated spin wave 
frequencies are a little stiffer than those found experimentally. 

As remarked above, in Fig. H] we show only the acoustical spin wave mode frequency 
regime. In Fig. El for the spectral densities in the innermost layer (upper panel) and the 
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FIG. 4: (color online) The spectral density in the innermost layer, in the acoustic spin wave regime, 
for wave vectors of Q\\ = +0.5A 1 (black curve) and Qy = — 0.5A 1 (red curve). Model C has 
been used for the calculation. In the ground state, the magnetization lies in plane along the long 
axis. 

outermost layer (lower panel) we show the spectral densities for the entire spin wave regime, 
including the region where the optical spin wave is found. It is clear that the spin orbit 
induced frequency shifts are largest for the optical mode which, unfortunately is not observed 
in the experiments.^ 

In Fig.[6]for a sequence of wave vectors, all chosen positive, we show a sequence of spectra 
calculated for the entire frequency range so both the acoustic and optical spin wave feature 
are displayed. The black curves show the spectral density of the innermost Fe layer, and 
the red curves are for the outer layer. The optical spin wave mode, not evident in the data, 
shows clearly in these figures. Notice that for wave vectors greater than lA 1 the acoustical 
mode is localized in the outer layer and the optical mode is localized on the inner layer. 
The optical mode is very much broader than the acoustical mode at large wave vectors, by 
virtue of the strong coupling to the electron hole pairs in the W 5d bands. 

An interesting issue is the absence of the optical mode from the SPEELS spectra reported 
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FIG. 5: (color online) For the wave vector Qu = 0.5A 1 we show the spectral densities in the 
innermost Fe layer (upper panel) and in the outermost layer (lower panel) for the Fe bilayer on 
W(110). The figure includes the optical spin wave feature. As in Fig.[H the black curve is calculated 
for Qu positive, and the red curves are for Qn negative. The calculations employ model C. 



in refs. |22| and |39|. We note that these spectra are taken with only two beam energies, 4 eV 
and 6.75 eV. At such very low energies, the beam electron will sample both Fe layers, 
so the SPEELS signal will be a coherent superposition of electron waves backscattered 
from each layer; the excitation process involves coherent excitation of both layers by the 
incident electron. As a consequence of the 180 degree phase difference in spin motions 
associated with the two modes it is quite possible, indeed even probable, that for energies 
where the acoustical mode is strong the intensity of the optical mode is weak, by virtue 
of quantum interference effects in the excitation scattering amplitude. In earlier studies 
of surface phonons, it is well documented that on surfaces where two surface phonons of 
different polarization exist for the same wave vector, one can be silent and one active in 
electron loss spectroscopy.- It would require a full multiple scattering analysis of the spin 
wave excitation process to explore this theoretically. While earlier— calculations that address 
SPEELS excitation of spin waves described by the Heisenberg model could be adapted for 
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Q (meV) 

FIG. 6: (color online) For the Fe bilayer and for several values of the wave vector (all positive), we 
show the spectral density functions for the innermost layer adjacent to the substrate (black curve) 
and those for the outer layer of the film. The calculations employ model C. 
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FIG. 7: (color online) For the Fe bilayer with magnetization in plane, we show the spin wave 
dispersion curves calculated with spin orbit coupling (black points) and without spin orbit coupling 
(red curves). Model C has been employed for these calculations. 

this purpose, in principle, a problem is that at such low beam energies it is necessary to take 
due account of image potential effects to obtain meaningful results.— This is very difficult 
to do without considerable information on the electron reflectivity of the surface.— It would 
be of great interest to see experimental SPEELS studies of the Fe bilayer with a wider range 
of beam energies to search for the optical mode, if this were possible. 

In Fig. [7J we show dispersion curves for the optical and acoustic spin wave branches for 
the bilayer. The magnetization lies in plane, and one can see that on the scale of this figure, 
the spin orbit effects on the dispersion curve are rather modest compared to those in the 
monolayer. For small wave vectors, with spin orbit coupling present, the dispersion curve of 
the acoustic spin wave branch is fitted by the form Q(Q\\) = 0.49 — 0.85Q\\ + 243Qy(meV) 
so at long wavelengths the influence of the Dzyaloshinskii Moriya interaction is more than 
one order of magnitude smaller than it is in the monolayer. 

If the magnetization is perpendicular to the surface, then symmetry considerations show 
that there are no left/right asymmetries in the spin wave propagation characteristics. One 
may see this as follows. Consider a wave vector Q\\ in the plane of the surface, which 
also is perpendicular to the magnetization, and thus perpendicular to the long axis. The 
reflection R in the plane perpendicular to the surface and which contains the magnetization 
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FIG. 8: (color online) For the bilayer and the case where the magnetization is perpendicular to the 
surface (model B), and for Qu = 0.6A , we show spectral density function calculated for positive 
values of Qu (continuous lines) and negative values of Qu (symbols). The black curve is the spectral 
density for the inner layer, and the red curve is the outermost Fe layer. 

simultaneously changes the sign of wave vector and the magnetization. If this is followed 
by the time reversal operation T, then Qu remains reversed in sign but the magnetization 
changes back to its original orientation. Thus the product RT leaves the system invariant 
but transforms Qu into —Q\\- The two propagation directions are then equivalent. 

We illustrate this in Fig. [8] where, for Qu = 0.6A , where it is shown that the spectral 
densities calculated for the two directions of propagation are identical, with spin orbit cou- 
pling switched on. Model B, in which the magnetization is perpendicular to the surface, has 
been used in these calculations. The spectral densities calculated for the two signs of Q\\ 
cannot be distinguished to within the numerical precision we use. 



IV. CONCLUDING REMARKS 



We have developed the formalism which allows one to include the influence of spin orbit 
coupling on the spin excitations of ultrathin ferromagnets on semi infinite metallic substrates. 
Our approach allows us to calculate the full dynamic susceptibility of the system, so as 
illustrated by the calculations presented in section III we can examine the influence of spin 
orbit coupling on the linewidth (or lifetime) of spin excitations, along with their oscillator 
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strength. As in previous work, we can then construct effective dispersion curves by following 
peaks in the spectral density as a function of wave vector, without resort to calculations of 



large numbers of very small distant neighbor exchange interactions. 
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'he results presented 
as discussed above, 



in Fig. H] are very similar to the experimental data reported in ref. 
though we see that in the bilayer the influence of the Dzyaloshinskii Moriya interaction is 
considerably more modest than in the monolayer. 

We will be exploring other issues in the near future. One interest in our minds is the 
influence of spin orbit coupling on the spin pumping contributions to the ferromagnetic 
resonance linewidth, as observed in ferromagnetic resonance (FMR) studies of ultrathin 
films. 43 It has been shown earlier 44 that the methodology employed in the present paper 
(without spin orbit coupling included) can be applied to the description of the spin pumping 
contribution to the FMR linewidth, and in fact an excellent quantitative account of the 
data on the Fe/Au(100) system was obtained. It is possible that for films grown on 4d 
and 5d substrates that spin orbit coupling can influence the spin pumping relaxation rate 
substantially. This will require calculations directed toward much thicker films than explored 
here. The formalism we have developed and described in the present paper will allow such 
studies in the future. 
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Appendix 

In this Appendix we provide explicit expressions for the various quantities which enter 
the equations displayed in Section II. While these expressions are unfortunately lengthy, it 
will be useful for them to be given in full. 
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A^ vl (U';mmf) = (A.1) 

The various expectation values in the equations above and those displayed below are 
calculated from the single particle Greens functions once the self consistent ground state 
parameters are determined. Then 

Bll tll > v >(ll';mm') = (Uwwv'icl^Ci'v^SirnSim' - Ui>^ u , vu >(cl^ci> v ^}5i> m 5ii m ^ 
Bll jlx ' u '(ll';mm') = Y (Ui'w'vyv( c lrt c i'nl)fii'rnfii'm' — Ui- Vfl '^ u '(cj v1: ci' u i}5i m 5i m i + 

V 

^11^^(11'; mm') = Y {Ui'.^^y^cj^Ci'^Si'rnSiirn' - U t . m > ^ v > (c\^C V v ±) 5i m 8i m > - 

—Ui>y U)VU >(c}^ci> v {)5[i m 5[i m ^ 
B% ill , vl (ll'-,mrri) = 

(A.2) 



V 

Bf u y u i{U'] mm 1 ) = Y [(Ui'w'vyn ~ U v -^ v ^ v i) {c\^c v ^)5 Vm 5 Vm/ - 
v 

— (Ui- W '^ U ' — Uiy^yi) {c lr) ^ci' V ^)8i m 8i m i 

Bf u y v '{U']mm') = Y {Uiy v y v (c\^c Vv ^)5 Vm 5 Vm i - ^ ;w ' iM ^(c]^qv t )^ m (^ m ') 
v 

B 2 ^y v i{U']mm') = -YUiiy u>vu i(c\^c Uri ])5v m 5v m > 

v 

(A.3) 



Bf u y u/ {U';mm') = -^^^/(c^Q/^cWW 

v 

B^ v y u i{ll']mm') = Y (Ui- tl i u y v (cl fli ci' v i}5i> m Si> m > - ^; W >^(c^Qv;)£z m (W) 

V 
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B™ tlJ i v i{ll'',rnm') - ^ \{U v .^ v y v - U v .^ v ^) {c\ pii c Vrii )5 Vm 5 Vm i- 
v 

Bll^ ul (ll';mm') =J2 U Wri,^'(4r r t C l'H) S lrn S lrn> 

(A.4) 



4v(«';-') = o 
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V 

( 



B]l^A ll 'i mm ') = T^5 lm 5^ - (T^ 1 ' t )*<W<W 
B^ >fl > v >(H';mm') = OL* v . vtv bi m bv m ib m i 



(A.5) 



Bli^AU';mm') = (A.6) 



B? u>ll/u ,(ll , ;mm , ) = 
Bli^ uf (U';mm') = -al^SimSum'S^ (A. 7) 
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B ltyu>( 11 '; mm> ) = tf'-yufilmfil'm'Sw' ( A - 8 ) 

mm') = frj5 lm 8^ - (T^)*<W<^ (A.9) 
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